// @HEADER
//*********************************************************************//
//  SiYuan: A numerical PDE solver                                     //
//  Copyright (2022) YUAN Xi                                           //
//  This Software is released under the BSD 2-Clause license detailed  //
//  in the file "LICENSE" in the top-level SiYuan directory            //
//*********************************************************************//
// @HEADER

#ifndef _BC_NEUMANN_IMPL_HPP
#define _BC_NEUMANN_IMPL_HPP

#include "Panzer_Workset.hpp"
#include "Panzer_GlobalEvaluationDataContainer.hpp"
#include "Phalanx_DataLayout_MDALayout.hpp"

// Evaluators
#include "Panzer_ConstantFlux.hpp"
#include "Panzer_Integrator_BasisTimesScalar.hpp"
#include "Panzer_Integrator_BasisTimesVector.hpp"
#include "Panzer_Constant.hpp"
#include "Panzer_ConstantVector.hpp"

#include "Teuchos_Array.hpp"
#include "Convection.hpp"
#include "Impedance.hpp"

namespace SiYuan {

namespace {

//　Fetch field value upon local cell basis from global ones 
template <typename ScalarT,typename LO,typename GO,typename NodeT>
class Assemble_Functor {
public:
  typedef typename PHX::Device execution_space;
  typedef PHX::MDField<const ScalarT,panzer::Cell,panzer::NODE> FieldType;

  Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> r_data;

  Kokkos::View<const LO**,PHX::Device> lids;    // local indices for unknowns
  Kokkos::View<const int*,PHX::Device> offsets; // how to get a particular field
  FieldType field;

  KOKKOS_INLINE_FUNCTION
  void operator()(const unsigned int cell) const
  {

    // loop over the basis functions (currently they are nodes)
    for(std::size_t basis=0; basis < offsets.extent(0); basis++) {
       int offset = offsets(basis);
       LO lid    = lids(cell,offset);
       Kokkos::atomic_add(&r_data(lid,0), field(cell,basis));

   } // end basis
  }
};

}


//**********************************************************************
template<typename EvalT, typename Traits>
NeumannEvalutor<EvalT, Traits> :: NeumannEvalutor( const NeumannBoundary& bc, 
  const Teuchos::RCP<const panzer::GlobalIndexer> & indexer,
  Teuchos::RCP<panzer::PhysicsBlock>& side_pb)
: PHX::EvaluatorWithBaseImpl<Traits>("Neumann Boundary Conditions"), 
  _need_scatter(false), globalIndexer_(indexer)
{
  m_sideset_name = bc.m_sideset_name;
  m_eblock_name = bc.m_eblock_name;
  //m_group_id = bc.m_group_id;
  m_dof_name = bc.m_dof_name;
//  m_bc_type = bc.m_bc_type;
  std::string val_type = bc.m_pl_value.get<std::string>("Type");
//  std::size_t found = val_type.find("Complex");
 // if (found!=std::string::npos)
 //   m_value = XYZLib::variable<std::complex<double>>(bc.m_pl_value);
 // else
    m_value = XYZLib::variable<double>(bc.m_pl_value);

//  const std::map<int,Teuchos::RCP< panzer::IntegrationRule > >& ir = side_pb.getIntegrationRules();
  int integration_order; // = side_pb.getIntegrationOrder();
  Teuchos::RCP<panzer::IntegrationRule> ir = this->buildIntegrationRule(integration_order, *side_pb);
  Teuchos::RCP<const panzer::PureBasis> pbasis = this->getBasis(m_dof_name, *side_pb);
  Teuchos::RCP<const panzer::FieldLayoutLibrary> fll = side_pb->getFieldLibrary()->buildFieldLayoutLibrary(*ir);
  Teuchos::RCP<panzer::BasisIRLayout> basis = fll->lookupLayout(m_dof_name);

  Teuchos::ParameterList p("Neumann BC");
  Teuchos::RCP< PHX::Evaluator<panzer::Traits> > op;
  try {
    if( m_dof_name == "Flux" ) {
    //  p.set("Name", m_flux_name);
      if( pbasis->isVectorBasis() )
      {
        p.set("Data Layout", ir->dl_vector);
     //   p.set("Value", params.get< Teuchos::Array<double> >("Value"));
     //   if(pbasis->dimension()>1)
     //     p.set("Value Y", this->m_bc.params()->template get<double>("Value Y"));
     //   if(pbasis->dimension()>2)
     //     p.set("Value Z", this->m_bc.params()->template get<double>("Value Z"));
    
        op = Teuchos::rcp(new panzer::ConstantVector<EvalT,panzer::Traits>(p));
      }
      else {
        p.set("Data Layout", ir->dl_scalar);
        p.set("Type", "Constant");
     //   p.set("Value",  params.get< double >("Value"));
        op = Teuchos::rcp(new panzer::Constant<EvalT,panzer::Traits>(p));
      }
    }
    else if( m_dof_name == "Convection" )
    {
    //  p.set("Name", m_flux_name);
    //  p.set("Ambient Temperature", this->m_bc.params()->template get<double>("Ambient Temperature"));
    //  p.set("Heat Transfer Coefficient", this->m_bc.params()->template get<double>("Heat Transfer Coefficient"));
      p.set("DOF Name", m_dof_name);
      p.set("IR", ir);
      p.set("Basis", basis); 
      op = Teuchos::rcp(new Convection<EvalT,panzer::Traits>(p));
    }
    else if( m_dof_name == "Radiate" )
    {
      p.set("Name", m_flux_name);
    //  p.set("Ambient Temperature", this->m_bc.params()->template get<double>("Ambient Temperature"));
    //  p.set("Radiate Coefficient", this->m_bc.params()->template get<double>("Radiate Coefficient"));
      p.set("DOF Name", m_dof_name);
      p.set("IR", ir);
      p.set("Basis", basis); 
      op = Teuchos::rcp(new Radiate<EvalT,panzer::Traits>(p));
    }
    else if( m_dof_name == "Acoustics Velocity" )
    {
    //  auto gd = this->getGlobalData();
    //  double omega = gd->pl->template getRealValue<EvalT>("Frequency");
    //  p.set("Name", m_flux_name);
    //  p.set("Resistance", this->m_bc.params()->template get<double>("Resistance"));
    //  p.set("Reactance", this->m_bc.params()->template get<double>("Reactance"));
    //  p.set("Density", this->m_bc.params()->template get<double>("Density"));
    //  p.set("Speed", this->m_bc.params()->template get<double>("Speed of Sound"));
      p.set("DOF Name", m_dof_name);
      p.set("IR", ir);
      p.set("Basis", basis); 
    //  p.set("Omega", omega);
      op = Teuchos::rcp(new AcousticVelocity<EvalT,panzer::Traits>(p));
    }
    else
      throw std::runtime_error("Neumann Strategy Not Defined!");
  }
  catch (std::exception& e) {
      std::cout << e.what() << std::endl;
  }

  m_residual_name = "Residual_" + m_dof_name;

  Teuchos::RCP<PHX::DataLayout> dummy = Teuchos::rcp(new PHX::MDALayout<void>(0));
  const PHX::Tag<EvalT> fieldTag("Neumann Condition", dummy);
  this->addEvaluatedField(fieldTag);
  this->setName(m_residual_name);
};


template<typename EvalT, typename Traits>
void NeumannEvalutor<EvalT, Traits> :: preEvaluate(typename Traits::PreEvalData d)
{
  if(Teuchos::is_null(m_GhostedContainer))
    m_GhostedContainer = Teuchos::rcp_dynamic_cast<panzer::LinearObjContainer>(d.gedc->getDataObject("Ghosted Container"));

//  if(m_GhostedContainer==Teuchos::null) {
    // extract linear object container
//    m_GhostedContainer = Teuchos::rcp_dynamic_cast<panzer::LOCPair_GlobalEvaluationData>
//              (d.gedc->getDataObject("Ghosted Container"),true)->getGhostedLOC();
//  }
}

// ***********************************************************************
/*template <typename EvalT>
void NeumannEvalutor<EvalT>::
postRegistrationSetup(typename panzer::Traits::SetupData d, PHX::FieldManager<panzer::Traits>& fm)
{
  if( !this->isActive() ) return;
}
*/
// ***********************************************************************
template<typename EvalT, typename Traits>
Teuchos::RCP<panzer::PureBasis> NeumannEvalutor<EvalT,Traits>::
getBasis(const std::string dof_name,const panzer::PhysicsBlock& side_pb) const
{
  const std::vector<std::pair<std::string,Teuchos::RCP<panzer::PureBasis> > >& dofBasisPair = side_pb.getProvidedDOFs();
  Teuchos::RCP<panzer::PureBasis> basis;
  for (std::vector<std::pair<std::string,Teuchos::RCP<panzer::PureBasis> > >::const_iterator it = 
	 dofBasisPair.begin(); it != dofBasisPair.end(); ++it) {
    if (it->first == dof_name)
      basis = it->second;
  }
  
  TEUCHOS_TEST_FOR_EXCEPTION(is_null(basis), std::runtime_error,
			     "Error the name \"" << dof_name
			     << "\" is not a valid DOF for the boundary condition:\n"
			     << this->m_bc_type << "\n");

  return basis;
}

// ***********************************************************************
template<typename EvalT, typename Traits>
Teuchos::RCP<panzer::IntegrationRule> NeumannEvalutor<EvalT,Traits>::
buildIntegrationRule(const int integration_order,const panzer::PhysicsBlock& side_pb) const
{
  TEUCHOS_ASSERT(side_pb.cellData().isSide());
  Teuchos::RCP<panzer::IntegrationRule> ir = Teuchos::rcp(new panzer::IntegrationRule(integration_order,side_pb.cellData()));
  return ir;
}


template<typename EvalT, typename Traits>
void NeumannEvalutor<EvalT, Traits> :: evaluateFields(typename Traits::EvalData d)
{
  if(_need_scatter)
  {
   /* auto x_data = m_GhostedContainer->getLocalView_x<PHX::Device>();

    globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);

    auto lids = scratch_lids_;
    for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
      auto offsets = scratch_offsets_[fieldIndex];
      auto gather_field = gatherFields_[fieldIndex];

      Kokkos::parallel_for(localCellIds.size(), KOKKOS_LAMBDA (std::size_t worksetCellIndex) {
       // loop over basis functions and fill the fields
        for(std::size_t basis=0;basis<offsets.extent(0);basis++) {
          int offset = offsets(basis);
          LO lid    = lids(worksetCellIndex,offset);

          // set the value and seed the FAD object
          gather_field(worksetCellIndex,basis) = x_data(lid,0);
        }
      });
    }*/
  }
  /* std::vector<panzer::LocalOrdinal> indx;
  for (auto i = m_Neumanns.begin(); i != m_Neumanns.end(); ++i){
    indx.emplace_back(i->first);
  }
  m_GhostedContainer->evalNeumannResidual(m_Neumanns); */
}


//**********************************************************************
template<typename EvalT, typename Traits>
CLoadEvalutor<EvalT, Traits> :: CLoadEvalutor( const std::vector< CLoad >& bcs, 
  const Teuchos::RCP<const panzer::GlobalIndexer> & indexer )
: PHX::EvaluatorWithBaseImpl<Traits>("CLoad Boundary Conditions"), 
   globalIndexer_(indexer), called(0)
{
  m_cloads.clear();
  for( auto bc: bcs )
  {
  //  if( !bc.isActive() ) continue;
    double v = bc.m_params.get< double >("Value");
  //  m_nodeset_name = bc.m_nodeset_name;
    //m_group_id = bc.m_group_id;
    m_dof_name = bc.m_dof_name;

    const int fnum = globalIndexer_->getFieldNum(m_dof_name);
    globalIndexer_->getNodesetsLocalIndex( fnum, bc.m_nodes_gid, m_ldofs);

    for( auto dof: m_ldofs )
    {
      m_cloads.insert( std::make_pair(dof, v) );
    }
  }

  Teuchos::RCP<PHX::DataLayout> dummy = Teuchos::rcp(new PHX::MDALayout<void>(0));
  const PHX::Tag<double> fieldTag("Cload Condition", dummy);
  this->addEvaluatedField(fieldTag);
  this->setName("Cload_Residual");
};

template<typename EvalT, typename Traits>
void CLoadEvalutor<EvalT, Traits> :: preEvaluate(typename Traits::PreEvalData d)
{
  if(Teuchos::is_null(m_GhostedContainer))
    m_GhostedContainer = Teuchos::rcp_dynamic_cast<panzer::LinearObjContainer>(d.gedc->getDataObject("Ghosted Container"));
}

template<typename EvalT, typename Traits>
void CLoadEvalutor<EvalT, Traits> :: evaluateFields(typename Traits::EvalData d)
{
//  std::cout << "RES:" << PHX::print<EvalT>()  << "," << called << std::endl;
  if( called<2 ) {
    m_GhostedContainer->applyConcentratedLoad(m_cloads);
    called++;
  }
}

}

#endif
